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Abstract 

In  this  paper,  numerical  methods  are  proposed  for  Poisson  equations  defined  in  a  finite  or  in 
the  infinite  domain  in  three  dimensions.  In  the  domain,  there  can  exists  an  interface  across  which 
the  flux  and  the  solution  are  discontinuous.  To  deal  with  the  discontinuity  in  the  source  term 
and  in  the  flux,  the  original  problem  is  transformed  to  a  new  one  with  a  smooth  solution.  Such 
a  transformation  can  be  carried  out  easily  through  an  extension  of  the  jumps  along  the  normal 
direction  if  the  interface  is  expressed  as  the  zero  level  set  of  a  three  dimensional  function.  An 
auxiliary  sphere  is  used  to  separate  the  infinite  region  into  an  interior  and  exterior  domain.  The 
Kelvin’s  inversion  is  used  to  map  the  exterior  domain  into  an  interior  domain.  The  two  Poisson 
equations  defined  in  the  interior  and  the  exterior  written  in  spherical  coordinates  are  solved 
simultaneously.  By  choosing  the  mesh  size  carefully  and  exploiting  the  fast  Fourier  transform, 
the  resulting  finite  difference  equations  can  be  solved  efficiently.  The  approach  in  dealing  with 
the  interface  has  also  been  used  with  the  artificial  boundary  condition  technique  which  truncates 
the  infinite  domain.  Numerical  results  demonstrate  second  order  accuracy  of  our  algorithms. 

Key  words,  arbitrary  interface,  fast  3D  Poisson  solver,  immersed  interface  method,  infinite 
domain,  extension  of  jumps,  spherical  coordinates,  level  set  function,  artificial  boundary  condition. 
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1  Introduction 


In  this  paper,  we  consider  the  three-dimensional  Poisson  equation  of  the  form 


A«(x) 


V-Mi(x),  if  (x)  eft", 

V  •  M2(x),  if  (x)  <E  =  M3  -Q-, 


lim  w(x)  =  0. 

|x|— KX> 


(1.1) 


We  assume  that  Mj(x),  1  <  i  <  2,  is  continuous  and  bounded  in  each  domain  of  its  definition.  In 
many  applications,  we  simply  have  M2(x)  =  0.  Across  the  interface  between  $2“  and  Q+,  there  is 
a  finite  jump  in  the  flux  of  the  solution 


'U'n  — 


du 

dn 


—  [  Vu  ■  n]  =  [  (M2  —  Mi)  •  n] 


(1.2) 
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where  n  is  the  unit  normal  direction  pointing  outward  along  the  interface  T  —  .  The  jump  is 

defined  as  the  difference  of  the  limiting  values  of  the  quantity  between  the  fl+  side  and  the  fl-  side. 
Figure  1  illustrates  the  description  of  the  problem  in  the  case  of  M2  =  0  and  V  •  Mi  =  /(x).  The 
jump  conditions  can  be  derived  either  through  physical  reasoning  or  mathematical  derivations,  see 
[5,  11,  12]. 


A  u  —  0 


Figure  1:  A  diagram  of  the  problem  with  M2  =  0  and  V  •  Mi  =  /(x)  in  a  bounded  domain  fl_. 
The  problem  is  defined  in  the  entire  space.  An  auxiliary  sphere  r  —  a  is  used  to  divide  the  infinite 
domain  into  two  parts. 

Modeling  and  simulation  of  ferromagnetic  materials,  see  [5,  3],  is  one  of  many  applications  of 
the  problem  stated  above.  Ferromagnetic  materials  are  widely  used  as  recording  media.  They  are 
also  currently  being  explored  as  an  alternative  to  semi-conductors  as  memory  devices.  Let  M  be 
the  magnetization  vector  field,  then  the  magnetic  field  is  given  by  —  Vu,  where  the  potential  u 
satisfies  the  Poisson  equation 

A«  =  V-(Xxn-),  (1-3) 

where  Xn~  is  fhe  characteristic  function  of  fl-.  Here,  fl-  is  the  domain  occupied  by  the  ferromag¬ 
netic  material,  see  [5]  for  a  detailed  description.  Therefore,  the  source  term  has  a  delta  function 
singularity  which  causes  a  discontinuity  in  the  flux  of  the  solution  across  the  boundary  of  fl-.  The 
domain  fl-  is  typically  a  closed  surface  (not  necessarily  a  sphere). 

It  is  important  to  provide  efficient  and  reliable  numerical  solutions  for  the  equation  above. 
As  quoted  from  [5],  “Even  though  there  is  an  extensive  literature  on  numerical  simulations  using 
micro -magnetics,  see  [1]  for  a  review,  with  even  commercially  available  software,  the  accuracy  of 
these  numerical  computations  is  still  a  serious  issue.” 

In  this  paper,  we  propose  second  order  accurate  fast  algorithms  based  on  the  fast  Fourier 
transform  (FFT)  and  the  immersed  interface  method  (IIM)  [10,  12,  15,  13,  16]  to  solve  the  three 
dimensional  Poisson  equation  with  discontinuities  in  the  flux  and  the  source  term  in  the  entire 
space. 

To  solve  (1. 1)- (1.2)  numerically,  we  need  to  overcome  two  difficulties  associated  with  the  prob¬ 
lem.  The  first  one  is  how  to  solve  the  Poisson  equation  in  the  infinite  domain.  The  second  one  is 
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how  to  handle  the  discontinuities  in  the  source  term  and  in  the  flux  of  the  solution. 

Finite  difference  approaches  are  often  preferred  for  a  number  of  reasons.  They  are  simple  to 
implement  and  fast  when  coupled  with  FFT.  For  non-linear  problems  in  which  a  Green  function 
may  not  be  available,  or  for  free  boundary  /moving  interface  problems,  finite  difference  methods 
with  fixed  underlying  grids  are  more  attractive.  The  challenge  is  how  to  maintain  the  accuracy 
near  the  boundary  or  the  interface. 

For  non-interface  problems  defined  in  an  infinite  domain,  we  refer  the  readers  to  [20]  for  a 
review  on  the  regularity  of  the  solution  and  a  variety  of  numerical  methods,  particularly,  the 
artificial  boundary  condition  (ABC)  techniques.  Using  ABC  techniques,  one  needs  to  truncate  the 
infinite  domain  and  then  to  impose  some  boundary  conditions  at  the  truncated  boundary.  In  [2], 
the  authors  derived  a  sequence  of  different  accurate  boundary  conditions  for  Helmholtz  equations 
along  the  artificial  sphere.  Some  other  artificial  boundary  conditions  can  be  found  in  [20]. 

In  this  paper,  we  propose  an  algorithm  that  can  solve  the  Poisson  equation  on  the  entire  three 
dimensional  space.  The  main  idea  is  to  use  an  auxiliary  sphere  to  separate  the  infinite  domain 
and  use  the  Kelvin’s  inversion  to  transform  the  exterior  domain  to  a  finite  interior  domain.  The 
mesh  points  corresponding  to  the  domain  outside  of  the  sphere  is  carefully  chosen  so  that  the 
finite  difference  discretizations  from  both  sides  can  match  together.  The  resulting  finite  difference 
equations  can  be  further  simplified  by  using  the  discrete  Fourier  transform  in  the  0  coordinate. 
Then  the  generalized  cyclic  reduction  method  can  be  used  to  solve  the  remaining  block  tri-diagonal 
linear  system.  Comparing  with  the  ABC  approach  that  solves  the  solution  only  inside  a  sphere, 
our  new  method  can  get  an  approximate  solution  in  the  entire  space,  but  the  computation  cost  is 
not  necessarily  doubled.  More  importantly,  we  do  not  have  to  worry  about  how  big  the  artificial 
sphere  is  or  what  kind  of  boundary  conditions  should  be,  which  makes  the  algorithm  potentially 
more  efficient. 

To  deal  with  the  discontinuity  across  the  interface  T  —  dil~ .  we  generalize  the  method  developed 
in  [16]  for  two  dimensional  problems  to  spherical  coordinates.  The  main  idea  is  to  extend  the  jumps 
in  the  solution  and  the  flux  along  the  normal  line  of  the  interface  in  a  neighborhood  of  T.  Then, 
we  can  construct  a  new  function  that  has  the  same  jump  conditions  as  the  solution  of  (1.1)- 
(1.2).  Using  a  transformation,  we  get  a  new  interface  problem  whose  solution  is  smooth  across  the 
interface.  Therefore,  a  second  order  accurate  scheme  can  be  obtained  using  the  standard  central 
finite  difference  scheme  in  spherical  coordinates  with  minor  modifications.  The  key  to  the  success 
of  this  approach  is  how  to  get  the  orthogonal  projections  of  the  grid  points  (in  the  neighborhood 
of  the  interface)  on  the  interface  which  is  not  trivial  when  the  geometry  is  described  in  spherical 
coordinates. 

Other  possible  approaches  include  the  use  of  an  integral  equation  like  the  fast  multipole  method 
[6] ,  since  the  solution  to  the  Poisson  equation  can  be  expressed  as  a  convolution  of  the  fundamental 
solution  with  the  source  term  /.  For  the  ferromagnetic  materials  simulation,  a  boundary  correction 
method  based  on  the  volume  integral  has  been  used  in  [3]. 

The  paper  is  organized  as  follows.  In  Sec.  2,  we  describe  the  transformation  based  on  an 
orthogonal  extension  of  known  jumps.  The  solution  to  the  new  Poisson  problem  is  smooth.  An 
efficient  fast  Poisson  solver  in  the  whole  space  and  the  artificial  boundary  condition  technique  are 
explained  in  Sec.  3.  The  entire  algorithm  for  the  interface  problem  defined  in  a  finite  or  the  infinite 
domain  is  outlined  in  Sec.  4.  Numerical  examples  and  grid  refinement  analysis  are  provided  in 
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Sec.  5. 


2  A  transformation  based  on  the  orthogonal  extension  of  the  jumps 


Without  loss  of  generality,  we  consider  the  following  more  general  problem 

fi  (x) ,  if  xefoj 
/2(x),  if  x  G  =  M3  —  fl_, 


Aa(x)  = 


u 


w , 


\Ur, 


lim  u(x)  =  0 

jx|— >oo 


(2.4) 


where  w  G  C2(r)  and  u  G  C2{T)  are  two  functions  defined  only  on  the  interface  T  —  dil~ .  First, 
we  briefly  review  how  to  transformation  the  above  interface  problem  to  a  problem  with  a  smooth 
solution.  This  was  first  introduced  in  [16]  for  two  dimensional  problems. 

Let  <p(x)  be  a  real-valued  function  such  that 


¥>(x) 


<  0,  if  x  G  fi~, 

—  0,  ifxGr  =  0fi",  (2.5) 

>  0,  if  x  G  fl+. 


We  assume  that  <p(x)  G  C3(M3 )  in  a  neighborhood  of  the  interface  T 1 .  so  that  the  zero  level  set 
of  ip,  <p(x)  —  0,  represents  the  interface.  Usually,  the  level  set  function  is  chosen  as  the  signed 
distance  function  (|V</j|  =  1)  from  the  interface,  see  [17,  18]  and  the  references  therein.  Therefore, 
it  is  reasonable  to  assume  that  |V<£>|  /  0  in  the  neighborhood  of  the  interface  T.  We  define  the 
extensions  of  w(X(s))  and  u(X(s))  along  the  normal  line  (both  directions)  as 

we(x)  —  we  (X(s)  +  an)  =  w(X(s)),  (2.6) 


and 


^e(x)  —  ve  (X(s)  T  an)  —  n(X(s)),  (2.7) 

for  all  a  G  R  such  that  the  normal  lines  do  not  intersect2,  where  n  is  the  unit  normal  direction 
pointing  outward.  We  construct  the  following  function  based  on  the  extension 

u(x)  =  we{x)  +  ve (x)  (2-8) 

Note  that  u(x)  G  C 2  in  the  neighborhood  of  the  interface  T  since  we  assume  that  "ic(s),  v(s)  are  in 

C2  in  the  domain  of  the  definition,  and  ip  is  in  C3  (M3 )  in  the  neighborhood  of  the  interface  T.  Let 

us  also  define 

r  o,  if  ip{x)  <  o, 

u{x)  —  H(ip(x))u(x)  —  <  |u(x),  if  <p(x)  —  0,  (2.9) 

_  1  ft(x),  if  tKx)  >  0, 

1In  implementation,  <p(x )  G  C2  seems  to  be  enough  for  second  order  accuracy. 

2 Theoretically,  there  is  always  a  neighborhood  that  the  normal  line  do  not  intersect.  Numerically,  in  case  the 
normal  lines  intersect,  we  can  simply  pick  up  an  extension  and  the  algorithm  still  works. 
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in  the  same  neighborhood  in  which  tt(x)  is  well  defined,  where  H(-)  is  the  Heaviside  function.  We 
have  the  following  theorem. 

Theorem  2.1  Letu(x)  be  the  solution  of  (2-4),  u(x)  be  defined  in  (2.9).  Define  q(x)  —  u(x)—u(x). 
Then  in  the  neighborhood  of  the  interface  where  we(x)  and  ve(x)  are  well  defined,  the  following  are 
true: 

Aq(x)  —  /(x)  —  H(tp(x))  Ait(x),  xrf-T,  (2.10) 

[q]  =  0,  [Vg-r]  =  0,  [Vg  •  n]  =  0,  (2.11) 

where  r  is  any  unit  tangent  direction  of  the  interface  T.  In  other  words,  the  new  function  q(x)  is 
a  sm.ooth  (C1)  function  across  the  interface  L\ 

The  proof  can  be  found  in  [16]. 

2.1  The  orthogonal  projection  on  the  interface  in  spherical  coordinates 

Using  Theorem  2.1,  we  can  transform  the  original  interface  problem  to  a  new  one  with  a  smooth 
solution.  In  order  to  take  advantage  of  the  transformation,  we  need  to  find  the  extensions  (2.6)-(2.7) 
along  the  normal  line  of  the  surface  T  to  get  u  and  u.  Since  our  fast  Poisson  solver  is  in  spherical 
coordinates,  it  is  natural  to  find  the  orthogonal  projection  of  a  point  near  the  interface  using  the 
information  (level  set  function  and  the  underlying  grid)  in  spherical  coordinates. 

The  spherical  coordinates  system  is 

x  —  r  sin  <f>  cos  9 , 

y  —  r  sincf  sin#,  r  >  0,  0  <  8  <  2tt,  0  <  <f>  <  n. 

z  —  r  cos  (f). 

Let  x  be  a  point  near  the  interface  T,  and  x*  be  the  corresponding  orthogonal  projection  of  x  on 
the  interface,  both  in  spherical  coordinates.  We  can  write 

x*  =  x  +  ap,  (2.12) 

where  a  is  a  scalar,  and  p  is  a  direction  to  be  determined.  If  the  direction  p  is  known,  we  can  use 
the  following  quadratic  equation 

<p(x)  +  (V<p(x)  ■  p)  a  +  ^  (p THe(<p(x))  p)  a2  =  0,  (2.13) 

to  approximate  the  scalar  a,  where  Vtp,  and  the  Hessian  matrix  He(<p)  are  evaluated  at  x.  The 
above  equation  is  a  third  order  accurate  approximation  to  <^(x*)  =  0.  The  key  here  is  how  to  decide 
the  right  direction  p  in  order  to  get  an  approximation  of  the  orthogonal  projection  x*. 

Let  [Aa;,  A y,  A z  ]T  and  [Ar,  Afi,  A 9  ]T  be  the  first  order  approximation  of  the  vector  x  —  x*  in 
Cartesian  and  spherical  coordinates,  respectively.  Then  we  have 


Aa: 

sin  <f>  cos  9  r  cos  cos  9  —r  sin  (f  sin  9 

Ar 

Ay 

sin  <f>  sin  9  r  cos  f>  sin  9  r  sin  cos  9 

Afi 

A  z 

cos  <f  —r  sin  0 

A  9 
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The  columns  of  the  Jacobian  matrix  above  are  orthogonal  and  the  Euclidean  norms  are  1,  r,  and 
|rsin^|,  respectively.  Therefore,  we  get 

d2  —  \\(x,y,z)  -  (x*  ,y* ,  z*)\\2  «  (Aa:)2  +  (Ay)2  +  (A  z]2  —  (A  r)2  +  (r  A^)2  +  (rsin</>A0)2. 

Since  <p(r*,  4>*,9*)  —  0,  we  have  p(r,  </>,  6)  —  p(r,  <^>,  8)  —  <p(r* ,  cf>* ,9*)  and 

v{rA,9)  -  (pr  Ar  +  (p^  A(f)  +  (pg  A8 

—  ipr  Ar  +  —  (r  A<fi) -\ - -  (r  sin  cj)  A8) ,  ^  ^ 

r  r  smip 

where  we  have  dropped  the  higher  order  terms.  Using  the  Cauchy-Schwartz’s  inequality,  we  have 

M  <  ]J{Vr)2  +  (^7)  +  '  V ^ Ar^2  +  ^  A^2  +  ^  sin  ^  A0^2' 

or  in  the  limit  case,  we  have 

j^-j-  <  \J (Ar)2  +  ( r  Ac/))2  +  ( r  sin  </>  A#)2. 

The  minimum  of  d  can  be  reached  if  the  equal  sign  can  be  reached  since  |<^|/V(/j|  is  an  approximation 
to  the  distance  between  (r,  4>1 9)  and  its  orthogonal  projection  (r* ,  (\*.  9*).  The  equal  sign  can  be 
reached  if  and  only  if  the  two  vectors  in  the  inner  product  (2.14)  are  collinear,  i.e.,  there  is  a  scalar 
a  such  that 


r  —  r  —  a  pr, 


,  ,*  V<t> 

<t>~  $  = 

9-9*  =  a 


Vo 


(r  sin  <fi) 2 

Thus  we  concluded  that  the  projection  is  along  the  direction 

-|T 


V<j> 

Vr,  9  5 


VO 


r 2  ’  (r  sin  (f))'2  \  ’ 

with  an  undetermined  scalar  factor  a  which  can  find  an  approximation  using  (2.13). 


(2.15) 


3  Fast  solvers  for  3D  Poisson  equation  in  the  infinite  domain 


The  Poisson  equation  in  the  three-dimensional  infinite  domain  It’’  can  be  written  in  spherical 
coordinates  as 


d2u  2  du  1  d2u  cot  du  1  d2u  _  . 

dr 2  r  dr  r2  d(j)2  r2  d(j)  +  r2  sin2  <j>  dff2  r’  ’ 


A  vanishing  condition  at  infinity 

u  — >  0  as  r  — >  00, 


(3.17) 
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is  imposed  in  order  to  ensure  the  uniqueness  of  the  solution. 

In  this  section,  we  will  introduce  two  different  fast  direct  methods  for  Eqs.  (3.16)-(3.17).  Both 
methods  involve  second-order  finite  difference  discretizations  to  the  Poisson  equation  (3.16)  which 
can  be  solved  efficiently  when  the  FFT-based  fast  solver  [9]  is  implemented.  The  only  difference 
of  the  two  methods  is  whether  the  problem  is  solved  in  the  whole  space  or  just  within  a  truncated 
domain. 

3.1  A  new  approach  without  truncation 

First,  we  present  a  new  approach  which  does  not  involve  any  truncation  of  the  infinite  domain.  Thus, 
a  numerical  solution  can  be  obtained  in  the  whole  three-dimensional  space.  Another  important 
advantage  of  this  approach  is  that  there  is  no  need  to  impose  any  artificial  boundary  conditions 
since  the  vanishing  condition  can  be  treated  naturally  in  our  finite  difference  setting.  Also,  the 
second-order  accuracy  is  preserved  in  the  whole  domain. 

In  [9],  Lai  et.  al.  have  developed  a  3D  FFT-based  fast  direct  solver  for  the  equation  (3.16) 
on  a  bounded  spherical  domain  {r  <  a}.  This  solver  has  three  major  features;  namely,  there  are 
no  pole  conditions,  a  fast  linear  algebraic  solver  can  be  applied,  and  the  scheme  is  able  to  handle 
different  boundary  conditions.  In  this  subsection,  we  will  develop  the  Poisson  solver  on  the  entire 
three-dimensional  space  based  on  the  method  on  a  finite  domain  described  in  [9]  and  the  Kevin’s 
transformation. 

Before  we  proceed,  let  us  decompose  the  infinite  domain  i?3  into  two  regions  using  a  sphere  which 
is  centered  at  the  origin  with  radius  a.  Let  the  interior  domain  be  —  {r  <  a}  and  the  exterior 
domain  be  =  {r  >  a},  respectively,  so  that  R 3  =  f2a  U  Then  the  boundary  d£la  will  be 
{r  —  a}.  The  present  idea  is  to  transform  the  part  of  Poisson  equation  in  the  exterior  unbounded 
domain  to  a  Poisson  equation  in  a  bounded  domain.  To  accomplish  that,  we  introduce  an 
inversion  mapping  f  =  a2/£,  or  (f,<j),6)  — >■  (a2/£,</>,  9)  that  maps  the  exterior  domain  {£  >  a}  into 
the  interior  domain  {f  <  a}.  Then,  we  define  the  functions  u  and  /  in  Qa  whose  values  are  related 
to  the  functions  u  and  /  in  by 

u{r,4>,0)  =  y  u  /(r,&0)  =  y  /  ,  f<a.  (3.18) 

The  above  transform  (3.18)  is  called  the  Kelvin’s  inversion.  The  functions  u  and  /  have  the  same 
regularities  in  the  domain  as  the  functions  u  and  /  in  iYa.  Also,  after  applying  the  Kelvin’s 
transform,  the  Laplace  operator  is  preserved  as 

2  1  cot  (f>  1 

+  JUC  +  -£2U4>4>  +  -|2 ~U4>  +  £2  sin 2^00 

f5(_  2  _  1  cot  4>  _  1 

—  n  l  Tt  T  ~Uf  T  T  zh  T  :  o  7^00 

a0  \  r  rz  rz  rz  sin  (.p 

One  can  directly  conclude  that  if  u  is  harmonic  in  iVa.  then  u  is  harmonic  for  every  points  in 
except  the  origin.  Therefore,  instead  of  solving  one  Poisson  equation  in  the  infinite  domain,  we 
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solve  the  following  two  Poisson  equations  in  the  finite  domain  Qa. 


2  1  cot  4>  1 

Urr  T  —Ur  T  ^ <b<b  T  o  “1“  — o  :  o  7  Uqq 
r  rz  rz  rz  sir  cp 

2  _  1  cot  <j)  _  1 

tiff  -)-  —  tif  +  To^<A<f>  +  —j  Ua  +  —  7  ^  T'U'dB 

r  rz  r2  r2  snr  <j> 


/,  in  fia, 

a4  f  •  0 

—  /,  m  na. 

r4 


(3.19) 

(3.20) 


The  above  two  equations  are  coupled  at  the  boundary  r  —  a  by  the  condition  u(a.  (j).  9)  —  a  u(a.  <j>,  9) . 
The  vanishing  boundary  condition  of  (3.17)  is  transformed  to  some  condition  of  u  at  the  origin, 
which  is  not  specified  explicitly.  We  shall  see  this  does  not  cause  any  trouble  in  our  finite  difference 
scheme  later.  Notice  that  the  values  of  u  and  u  and  their  derivatives  should  agree  with  each  other 
at  the  boundary  dUa  if  u  is  transferred  back  to  the  original  equation  and  the  geometry. 

We  now  discuss  how  to  solve  Eqs.  (3.19)-(3.20)  based  on  the  fast  direct  solver  developed  in  [9]. 
First,  we  choose  a  grid 


{ri,4>j,0k)  =  (*Ar,  (j  -  l/2)A<f>,kAdj  , 
where  the  mesh  spacings  for  each  direction  are 


1  <  i  <  M,  1  <  j  <  L,  1  <  k  <  N, 


(3.21) 


A  r  =  a/M,  A<j>  =  ir/L,  A  9  =  2t t/N. 


(3.22) 


In  order  to  match  the  position  of  the  ghost  grid  point  in  radial  direction,  the  mesh  width  Ar  is 
chosen  to  satisfy 

a2  ,  . 

a  —  Ar  — - — .  3.23 

a  +  Ar  v  ' 

This  requirement  has  the  meaning  that  the  ghost  grid  point  of  u  outside  the  boundary  {r  =  a}  will 

be  the  grid  point  of  u  which  is  the  closest  one  to  the  boundary.  By  a  simple  calculation,  we  obtain 

A f  =  a/(M  +  1).  So  the  grid  points  used  for  u  can  be  chosen  similarly  by 


(ri,(j)j,9k)  —  ^«Ar,  ( j  —  l/2)A</>, kA9^j  ,  1  <  i  <  M  +  1,  1  <  j  <  L,  1  <  k  <  N.  (3.24) 


Notice  that,  by  the  choice  of  Ar  and  Ar,  we  have  vm  —  tm+  i  =  a  and  thus  the  coupling  boundary 
condition  is  uM+i,j,k  =  auM,j,k- 

Once  we  have  set  up  the  grid,  we  can  apply  the  standard  second-order  centered  finite  difference 
approximation  to  the  equations  (3.19)  and  (3.20).  This  leads  to  the  following  finite  difference 
equations 


ui+l,j,k  ‘2‘uijk  +  ui—l,j,k  ,  2  'Wj— ,  1  ui,j+l,k  “^uijk  T  ui,j— l,k 

I  ^  A  I  O 


(Ar)2 


Vi 


2Ar  '  r2  (A<^>)2 


cot  (fij  'Ujjj+ljA;  ‘^’ijj  —  ljk  1 

I  O  rA  A  I 


2A  <j) 


r  2  sin2  (j)j 


'U‘i,j,k+ 1  “^Uijk  +  "Wj,  j,k—l 

(A  9)2 


=  fijk,  (3-25) 


ui+l,j,k  “^uijk  +  ui—l,j,k  ,  2  i,j,k  ui—l,j,k  ,  1  ui,j+l,k  ^uijk  T  l,fc 

“r  “  ,  .  “r  To 


(Ar)2 


Ti 


,  c°f  4*]  ui,j+l,k  ui,j—l,k  , 

>  —)  A  J  I 


2Af 

1 


f2 


(A^)2 


f2 


2A(f> 


r2  sin2  ^ 


“ijjfc+i  2 Uijk  +  Uijfi- 1  _  °4  /  /  9 

- (A0p - “  (3'26) 
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Since  the  the  solution  is  periodic  in  9  direction,  we  can  approximate  u  by  the  truncated  Fourier 
series 

N/  2-1 

u{ri,4>j,0k)  =  v,n(rii  fa)  ein9k }  (3.27) 

n=-N/  2 


where  un(r ,,  (f>3)  is  the  discrete  Fourier  coefficient  given  by 


N- 1 


u7 


{riAj)  =  ^ 


(3.28) 


k=0 


The  above  transformation  between  the  physical  space  and  Fourier  space  can  be  efficiently  performed 
using  the  FFT  with  0(N  log2  N)  arithmetic  operations.  The  expansions  for  the  functions  /,  u.  and 
/  are  similar. 

Substituting  those  expansions  into  the  standard  centered  finite  difference  scheme  (3.25)-(3.26) 
and  equating  the  corresponding  Fourier  coefficients,  the  finite  difference  equations  for  the  nth 
Fourier  coefficients  un  and  un  are  as  follows 


Ui+  2  Uij  +  Uf-u  2  Ui+i,j  Ui— ij  1  U,j-\ |-i  2  Uj  +  U,j—i 

(Ar)2  ri  2  A r  rf  (A</>)2 

cot  <j)j  U,j+ 1  -  U,j- 1  1  2  cos(nA0)  -2  _ 

r2  2A«£  r2  sin2  ^  (A0)2  ^ 


(3.29) 


U%-\- 2  C/jj  +  U—i,j  2  C/j+ij  Ui—ij  1  Uij+i  2  Uj  +  U,j— i 

(Ar)2  fj  2Ar  rf  (A</>)2 

cot  <j)j  Ujj+i  -  Ujj- 1  1  2  cos(tcA0)  -  2  -  _  =, 

rf  2A 4>  rf  sin2  cj)j  (A9)2  l]  rf 


(3.30) 


Here,  we  denote  the  discrete  values  by  Uij  —  un(ri,(f)j ),  Uj  —  un(ri,  (f)j).  Fij  —  fn(ri,<f>j)  and 
Fij  — 

We  now  discuss  how  to  handle  the  numerical  boundary  conditions  for  <j>  —  0  and  ({>  =  tt.  Using 
a  staggered  grid  in  the  <p  direction,  we  do  not  put  the  grid  points  on  the  north  (r/>  =  0)  and 
south  ((f)  —  7 r)  poles  directly.  Thus,  the  numerical  boundary  values  for  Uo  and  U,l+i  can  be 
easily  obtained  by  the  symmetry  constraint  of  Fourier  coefficients  [9];  that  is,  Uo  —  (— l)nUn  and 
U,l+ i  =  (-1  )nUiL-  Similarly,  we  have  U0  —  (~l)nUu  and  U^l+i  —  (-1  YUl-  Therefore,  no 
complicated  pole  treatments  are  needed  in  our  setting.  More  surprisingly,  the  numerical  boundary 
values  Uoj  and  Uoj  will  not  be  needed  either  since  the  coefficients  of  those  values  in  the  finite 
difference  equations  (3.29)-(3.30)  are  zero.  One  can  easily  check  this  from  the  above  finite  difference 
equations  in  the  case  of  i  —  1.  The  numerical  boundary  conditions  for  Umj  and  Um+ij  are  provided 
by  the  coupling  boundary  condition  Um+ij  —  a  Um3  ■  By  ordering  the  unknowns  appropriately, 
the  resulting  linear  system  of  equations  for  Uj  and  Uj  can  be  solved  altogether  by  the  generalized 
cyclic  reduction  algorithm  [19]  with  0(LM  log2M)  operations. 
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3.2  An  artificial  boundary  approach 


The  traditional  numerical  computation  for  (3.16)-(3.17)  is  to  truncate  the  unbounded  region  by  a 
sphere,  say  r  —  a;  then  an  artificial  boundary  condition  is  imposed  at  r  =  a.  The  resulting  problem 
is  then  discretized  and  solved  in  the  bounded  region  {r  <  a}.  However,  the  main  difficulty  is  how  to 
choose  the  artificial  boundary  condition  at  r  —  a  so  that  the  vanishing  condition  at  infinity  can  be 
translated  appropriately  to  the  boundary.  In  [2],  the  authors  developed  a  sequence  of  local  boundary 
conditions  for  some  elliptic  equations.  Those  artificial  boundary  conditions  can  be  considered  as 
applying  the  differential  operators  to  the  solution  at  the  boundary.  The  order  of  accuracy  of  those 
boundary  conditions  depends  on  the  order  of  the  highest  derivative  in  the  differential  operators.  In 
this  paper,  we  use  the  following  second-order  accurate  boundary  condition 


d2u  4  du  2u  _ 

Q__  T  2  ■ 


at  r  —  a. 


(3.31) 


dr2  '  r  dr 

The  detail  of  the  derivation  of  the  above  condition  can  be  found  in  [2].  In  short,  we  now  solve  the 
equation  (3.16)  in  a  finite  domain  {r  <  a}  with  the  boundary  condition  (3.31). 


The  finite  difference  discretization  for  the  resulting  problem  is  almost  the  same  as  the  one 
described  before.  That  is,  we  use  the  same  grid  as  described  in  (3.21)-(3.22)  and  follow  the  same 
solution  procedures  to  solve  one  Poisson  equation  in  the  domain  {r  <  a}.  The  treatments  for 
numerical  boundary  values  are  also  the  same  except  at  the  boundary  r  —  a  which  uses  the  condition 
(3.31)  instead.  Therefore,  the  numerical  boundary  value  of  Um+ i,j  can  be  found  by  the  second-order 
approximation  of  (3.31) 


Um+ i,j 


2  Um]  +  Um-i,]  4  Um+ i,j  ~  Um- i,j  2  XJuj  _  „ 

(Ar)2  rM  2  A r  rli 


(3.32) 


Once  we  incorporate  the  numerical  boundary  values  into  the  difference  scheme,  the  resulting  linear 
equations  are  again  solved  by  the  generalized  cyclic  reduction  method. 


Remark  3.1  If  the  domain  is  a  finite  r  <  a  instead  of  the  infinite  domain,  then  we  can  apply  the 
exact  boundary  condition  at  r  —  a  instead  of  the  artificial  boundary  condition. 

Let  us  close  this  section  by  summarizing  the  algorithms  and  the  operation  counts  in  the  following 
three  steps: 

1.  Compute  the  Fourier  coefficients  for  the  right-hand  side  function  as  in  (3.28)  by  FFT,  which 
requires  0(M LN  log2  N)  operations. 

2.  Solve  the  block  tridiagonal  linear  systems  for  N  Fourier  coefficients  by  the  generalized  cyclic 
reduction  algorithm,  which  requires  0{NLM\og2M )  operations. 

3.  Invert  the  Fourier  coefficients  as  in  (3.27)  by  FFT  to  obtain  the  solution,  which  requires 
0(MLN\og2N)  operations. 

The  overall  operation  count  is  roughly  0(NLM  log2(M  +  IV))  for  M  x  L  x  N  grid  points.  One 
should  notice  that  the  artificial  boundary  condition  approach  saves  about  the  half  of  the  work 
comparing  with  the  approach  without  truncation  since  the  former  does  not  compute  the  solution 
outside  the  truncating  boundary. 
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4  Dealing  with  the  interface  and  an  outline  of  the  algorithm. 


In  this  section,  we  explain  how  to  handle  the  interface  and  give  an  outline  of  our  algorithm.  For 
the  interface  problem  (2.4),  using  the  transformation  described  in  Section  2  and  a  level  set  function 
to  express  the  interface,  we  can  get  a  discrete  Poisson  equation  with  a  modified  right  hand  side. 
The  discrete  Poisson  equation  can  be  solved  either  in  the  entire  space  or  in  a  sphere  r  <  a  using 
an  artificial  boundary  condition  technique.  Below  we  explain  our  algorithm  step  by  step. 

Step  1:  Indexing  grid  points.  First  we  generate  a  grid  (3.21)-(3.22)  and  choose  a  sphere 
r  —  a  that  encloses  the  interface  T  —  dil~ .  A  level  set  function  whose  zero  level  set  (<^(r,  ip.  9)  —  0) 
is  the  interface  T.  is  defined  at  grid  points  as  fijk-  If  we  only  want  to  obtain  the  solution  in  a 
bounded  domain,  then  we  have  to  choose  the  artificial  sphere  r  —  a  carefully  in  order  to  use  an 
artificial  boundary  condition  technique  while  maintaining  second  order  accuracy. 


Let 


Vijk  —  max{^i— ijj.,  Vi,j,kt  Vi+l,j,ki  Vi ,j—  l,kt  Vi,j+ l,kt  Vi,j,k— 1:  Vi ,j,k+l\ i 
Vijk1  =  min-f^j-i^fc,  Vi,j,ki  Vi+l,j,ki  Vi,j-l,ki  Vi,j+ l,k:  Vi,j,k-li  Vi,j,k+ 1}- 

We  define  xt]k  as  an  irregular  grid  point  if 


(4.33) 


Vijk  '  Vijk  — 


min  ^  q 


(4.34) 


We  define  =  (r*,  (f)J:  9k)  as  an  sub-irregular  grid  point  if  it  is  not  an  irregular  grid  point,  but 
one  of  its  six  neighbors  is  an  irregular  grid  point.  If  Xjj^  is  neither  an  irregular  nor  a  sub-irregular 
grid  point,  then  we  call  it  a  regular  grid  point. 


Step  2:  Finding  the  orthogonal  projections.  If  Xijk  is  an  irregular  or  sub-irregular  grid 
point,  we  use  the  second  order  central  finite  difference  scheme  to  find  the  approximate  orthogonal 
direction  p  defined  in  (2.15)  at  Xjjyti  and  the  Hessian  matrix  He(ip(xt:jk))  in  spherical  coordinates. 
Then  we  solve  the  quadratic  equation  (2.13)  to  get  a  and  an  approximate  orthogonal  projection 
x*j{,  of  Xijk  on  the  interface  T. 

Step  3:  Extending  the  jumps.  We  need  to  extend  the  jumps  to  all  irregular  and  sub¬ 
irregular  grid  points  x^k  from  their  orthogonal  projection  x*.-k  on  the  interface  according  to  the 
following 

{we)ijk  -  w{x*ijk{s )),  (ve)ijk  -  v(x*jk{s)).  (4.35) 

Now  we  compute  Uijk  and  iiijk  defined  in  (2.8)  and  (2.9)  at  irregular  and  sub-irregular  grid  points 
Xijk  according  to  the  following 


fi ijk 


'U'ijk 


We(Xijk)+Ve(Xijk)  | v^(x^)|’ 

_  f  0,  if  v{xijk)  <  o, 

\  '9'ijki  if  v(xijk)  >  0, 


(4.36) 


Step  4:  Adding  the  correction  terms.  Since  the  solution  of  the  Poisson  equation  A (u(x)  — 
u(x)  —  f(x)  —  H(ip(x))  Aii(x)  has  a  smooth  solution,  an  approximate  solution  can  be  obtained  by 
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solving  the  following  finite  difference  equations 

fijk  ~  Hh{<fijk)Ahuijk  +  A huijk  +  Cijk,  if  xijk  is  irregular, 
fijk ,  otherwise. 


A h,Uijk  —  \ 


(4.37) 


Here,  C{jk  is  a  correction  term  needed  to  offset  the  discontinuities  in  second  order  derivatives,  and 
A/,  is  the  standard  discrete  Laplacian  in  spherical  coordinates  described  by 

a  TT  Ui+ijjk  2Uijk  +  Ui—ijjk  2  Ui- \-i,j,k  U{—  1  Ui,j+i,k  2Uijk  + 

—  - T~a  To - 1 - — - 1 - 


(Ar)2 

cot  4>j  Uj,j+ i,fc  ~  Uj,j-i,k 
r'2  2A  cj) 

and  Hk{<pijk)  is  the  discrete  Heaviside  function 


2  A  r 


rf  sin2  cf)j 


■i  m2 

Ui,j,k+1  —  2Uijk  +  k—l 

{Key2  : 


Hhi^Pijk) 


1  if  (pijk  >  0, 
0  if  Vijk  <  0. 


(4.38) 


The  correction  term  Cijk  is  defined  in  (4.40)  and  it  is  important  for  second  order  accuracy  of  our 
method.  Note  that  when  i  —  1,  the  coefficient  of  Uojk  in  the  expression  above  happens  to  be  zero. 
Therefore  we  have  avoided  the  singularity  at  r  —  0.  The  extension  of  the  jumps  to  sub-irregular 
grid  points  is  used  to  compute  AhUijk  and  Aku1]k  at  irregular  grid  points. 

Define 


^ijk  — 


The  correction  term  then  is 


fijk  ~  Hh(<pijk)Ahuijk ,  if  xijk  is  irregular, 
fijkl  otherwise. 


(4.39) 


Cijk—  Hh{  <Pi+ii,j+ji,k+ki  lPijk)'yi+ii,j+ji,k+kl 

k,iiM 


^i+ihj+Jhk+ki 

I  ^^i+ihj+jhk+k,  I 


JMuj+juk+ki 


~F, 


ijk ! 


(4.40) 


where  (ii,ji,kt)  =  {  (-1, 0, 0),  (1, 0, 0),  (0, -1,  0),  (0, 1, 0),  (0, 0, -1),  (0, 0, 1)  },  are  the 

coefficients  of  the  discrete  Laplacian  in  spherical  coordinates: 


7 'i±l,j,k  — 


(Ar)'2  r,Ar 
The  gradient  vector 


)  7»,j±l,fc  — 


V<£ 


1  ±  cot  4>j 

(ri  A 4>)2  2 r2  Acp 


i  li,j,k±l  — 


(ri  smcpjAQ)2 


Vri  7W  .  i 
r  r  sin  (p 


<Pe 


(4.41) 


(4.42) 


is  computed  using  the  standard  second  order  central  finite  difference  scheme  at  involved  grid  points. 

It  is  clear  now  that  the  discontinuity  in  the  solution  and  the  flux  only  affect  the  right  hand  side 
of  the  finite  difference  equations. 

Step  5:  Solving  the  discrete  system.  We  solve  the  discrete  system  of  equations  with 
modified  right  hand  side  using  the  fast  method  described  in  Section  3  to  get  an  approximate 
solution  to  the  original  problem  either  in  the  entire  space  or  inside  the  truncated  domain  r  <  a. 
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5  Performance  study 


We  have  done  intensive  tests  on  the  methods  discussed  in  this  paper  using  Sun  Ultra  workstations 
or  PCs.  Our  computer  codes3  have  not  been  optimized  and  parallelized.  However,  all  the  numerical 
results  confirm  second  order  accuracy  of  the  present  methods  in  the  infinity  norm. 


5.1  Example  1. 


In  this  example,  the  level  set  function  is 

ip(r,  <j>,9)  —  r2  —  (1  +  Br  sin  <j>  cos  9)  . 


(5.43) 


The  interface  in  Cartesian  grid,  therefore,  is  x2  +  y2  +  z2  —  Bx  —  1.  The  exact  solution  is  chosen  as 


The  source  term  is 


«(r,  4>,0) 


r3  sin  (j)  cos  6 
sin  <j>  cos  6 


if  tp(r,  4>.,  9)  <  0, 
if  <p(r,  4 f>,  9)  >  0. 


f{r,4>,9) 


lOrsin^cos#  if  ip(r,  (j),  9)  <  0, 
0  if  ip(r,  4>,  9)  >  0. 


(5.44) 


(5.45) 


The  jump  condition  is  determined  from  the  exact  solution  and  the  interface.  In  the  special  case 
when  B  —  0,  the  interface  is  the  unit  sphere  r  —  1,  the  jump  conditions  are 


[it]  =  0,  [un]  —  —  5  sin  ^  cos  0.  (5.46) 

If  B  /  0,  then  both  the  solution  and  the  flux  have  a  non-zero  jump. 

In  Table  1,  we  show  the  result  of  the  grid  refimenet  analysis  for  the  solution  both  inside  and 
outside  of  the  auxiliary  sphere  r  —  a  —  2.  The  error  is  defined  as 


E oo  =  max  {  | u(ri,  <pj,9k) 

ijk 


Uijk  |  }  •> 


(5.47) 


where  u(ri,(frj,9k)  and  U,.jk  are  the  exact  solution  and  computed  solutions  at  the  grid  point 
(ri.  (fij,  9k).  We  use  to  denote  the  maximum  error  inside  the  sphere,  and  E2^  to  denote  the 
maximum  error  outside  the  sphere.  The  order  of  accuracy  is  defined  as 


order  = 


\og(E(M)  /E(2M)) 
log  2 


(5.48) 


In  Table  1  (a),  B  —  0  and  the  interface  is  the  unit  sphere  which  is  aligned  with  the  grid  line 
in  r  direction.  The  solution  is  continuous,  but  both  the  flux  and  the  source  term  have  a  finite 
jump.  The  grid  refinement  analysis  show  clearly  second  order  accuracy  as  the  mesh  gets  finer.  In 
Table  1  (b),  the  interface  depends  on  r,  </>,  and  9.  It  is  not  aligned  with  any  grid  lines.  Now  the 
solution,  the  flux,  and  the  source  term  all  have  a  finite  jump  across  the  interface.  Generally,  for 
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Table  1:  The  grid  refinement  analysis  for  Example  1.  The  auxiliary  sphere  separating  the  infinite 
domain  is  r  —  a  —  2.  The  CPU  time  unit  is  in  seconds. 


(a)  Sphere  interface:  B  —  0.  The  solution  is  continuous  but  the  flux  is  not. 


M 

N 

L 

EL 

order 

EL 

order 

CPU(s) 

16 

32 

16 

7.0822  x  10“^ 

4.4105  x  10-3 

0.0900 

32 

64 

32 

1.7932  x  10-3 

1.9816 

1.1444  x  10-3 

1.9464 

0.6700 

64 

128 

64 

4.5196  x  10“4 

1.9883 

2.9192  x  10“4 

1.9710 

5.5200 

128 

256 

128 

1.1349  x  10“4 

1.9936 

7.3742  x  10“5 

1.9850 

49.870 

(b)  An  arbitrary  interface:  B  =  0.25.  Both  the  solution  and  the  flux  are  discontinuous. 


M 

N 

L 

EL 

order 

TP‘2 

•^oo 

order 

CPU(s) 

16 

32 

16 

5.0432  x  10-2 

9.0128  x  10-3 

0.0900 

32 

64 

32 

1.3027  x  10-2 

1.9529 

2.4347  x  10“3 

1.8882 

0.6600 

64 

128 

64 

1.8168  x  10_a 

2.8420 

6.3863  x  10“4 

1.9307 

5.5500 

128 

256 

128 

5.2148  x  10“4 

1.8008 

1.7160  x  10“4 

1.8960 

49.870 

Table  2:  The  grid  refinement  analysis  for  Example  1  using  the  artifical  boundary  condition  (3.31) 
and  B  —  0.25.  The  CPU  time  unit  is  in  seconds. 


M 

N 

L 

CM 

II 

Ls 

Kl 

order  ( a  —  2) 

EL(“  =  5) 

order  ( a  —  5) 

CPU(s) 

16 

32 

16 

BHiHlBI 

2.9658  x  10-2 

0.0700 

32 

64 

32 

2.4063 

6.8574  x  10-3 

2.1127 

0.4800 

64 

128 

64 

2.7957 

2.3333  x  10“3 

1.5553 

3.8200 

128 

256 

128 

2.0536 

6.2092  x  10“4 

1.9099 

32.070 
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interface  problems,  the  maximum  errors  usually  do  not  decrease  monotonically  (see  [14,  4]).  But 
the  average  convergence  rate  in  Table  1  (b)  shows  second  order  convergence  again. 

In  Table  2,  we  show  the  result  of  the  grid  refinement  analysis  using  the  artificial  boundary 
condition  (3.31).  The  error  is  measured  at  grid  points  inside  the  artificial  sphere  r  —  a.  The  CPU 
time  was  cut  by  less  than  half.  In  this  table,  we  list  the  numerical  results  for  a  —  2  and  a  —  5.  It  is 
hard  to  determine  the  optimal  value  of  a  since  the  choice  is  problem-dependent.  In  this  case,  the 
artificial  sphere  a  —  2  seems  to  be  a  good  choice  for  this  particular  artificial  boundary  condition 
technique.  When  we  increase  the  radius  of  the  sphere  to  a  —  5,  the  maximum  error  actually  gets 
a  little  bigger.  The  part  of  reason  is  that  the  mesh  resolution  when  a  —  5  is  worse  than  the  case 
when  a  —  2.  Nevertheless,  both  cases  show  second  order  accuracy. 

5.2  Example  2. 

In  this  example,  we  show  the  result  for  a  problem  where  the  solution  has  a  stronger  discontinuity. 
The  level  set  function  is 


(p(r,  <j>,  9)  —  r  —  (1  +  B  sin  (j>  cos  9) . 


The  exact  solution  is  chosen  as 


u(r,  <j>,  9) 


The  source  term  is 


rsin<^cos# 
2  sin  cf)  cos  9 


if  ip(r,  <j>,  9)  <  0, 
if  ip(r,  (j),9)  >  0. 


f{r,  4>,  0) 


— 15  r  sin  <^>cos  0  if  ip(r,  <j>,  9)  <  0, 
0  if  i p(r,  <j>,  9)  >  0. 


(5.49) 


(5.50) 


(5.51) 


The  jumps  in  [it]  and  \un]  are  determined  from  the  exact  solution  and  the  interface.  Note  that  even 
in  the  special  case  when  B  —  0,  the  solution  is  discontinuous.  In  the  following  tests,  we  choose 
B  =  0.1. 

In  Table  3,  we  show  the  grid  refinement  analysis  with  an  auxiliary  sphere  r  —  a  —  2  for  the 
solution  both  inside  and  outside  of  the  sphere.  Average  second  order  accuracy  is  again  verified. 

In  Table  4,  we  show  the  result  of  the  ABC  approach.  The  results  agree  with  our  analysis  for 
Example  1  and  Table  2. 

We  should  point  out  that  our  method  does  require  the  interface  to  be  smooth  enough  in  Eu¬ 
clidean  space  because  the  method  is  based  on  the  extension  of  the  jumps  along  the  orthogonal 
directions.  Some  interfaces  that  are  arbitrarily  differentiable  in  spherical  coordinates  may  have  a 
singularity  in  Euclidean  space.  For  example,  <p(r,  91  <f>)  —  r  —  (1  —  B  cos  9)  has  a  singularity  at  the 
pole  (j>  —  0  and  —  it  when  it  is  expressed  in  Cartesian  coordinates. 

3The  code  is  available  to  public  upon  request. 
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Table  3:  The  grid  refinement  analysis  for  Example  2.  The  auxiliary  sphere  separating  the  un¬ 
bounded  domain  is  r  =  a  —  2  and  B  —  0.1.  The  CPU  time  unit  is  in  seconds. 


M 

N 

L 

EL 

order 

-^oo 

order 

CPU(s) 

16 

32 

16 

IgsEEBHBil 

||l||g2|gj|sgj| 

0.080 

32 

64 

32 

0.8791 

3.0069 

0.6600 

64 

128 

64 

2.8282 

1.8998 

5.5500 

128 

256 

128 

1.6667 

1.9687 

50.090 

Table  4:  The  grid  refinement  analysis  for  Example  2  using  the  artifical  boundary  condition  (3.31). 
The  CPU  time  unit  is  in  seconds. 


M 

N 

L 

EL(a  =  2) 

order  (a  =  2) 

EL{a  =  5) 

order  (a  =  5) 

CPU(s) 

16 

32 

16 

KMEEBsarigi 

0.0700 

32 

64 

32 

2.2631 

3.4189 

0.4800 

64 

128 

64 

2.8602 

1.2003 

3.7800 

128 

256 

128 

1.6551 

2.0459 

31.9500 
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